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Abstract. The self-consistent interaction between a beam of charged particles and 
a wave is considered, within a Vlasov picture. The model is discussed with reference 
to the case of a Free Electron Laser. Starting with a spatially bunched waterbag 
distribution, we derive, via perturbative methods, closed analytical expressions for the 
time evolution of the main macroscopic observables. Predictions of the theory are 
shown to agree with direct numerical simulations. 
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1. Introduction 

Wave-particles interaction is a basic process in physics, which is encountered in a large 
class of different phenomena. Most spectacular examples are undoubtely found in 
astrophysical context, but also in laboratory plasmas where technological aspects of 
nuclear fusion are adressed. 

Free-Electron Lasers (FELs) [TJ [2] constitute another important application where 
the dynamical interplay between particles and waves is well known to be central. The 
physical mechanism that drives the process of light amplifiction is in fact the interaction 
between a relativistic electron beam and a co-propagating optical wave, inside the so- 
called undulator. The latter generates a magnetostatic periodic field, thus forcing the 
electrons to follow curved trajectories and emit synchrotron radiation. This incoherent 
light seeds, termed spontaneous emission, cumulate under resonance condition, and 
eventually result in the coherent laser signal. In a single-pass scheme, the laser is 
produced during a single passage inside the undulator, and the basic features of the 
system dynamics are successfully captured by a simple one- dimensional Hamiltonian 
model [5]. 

Remarkably, an analogous mathematical formulation is shown to describe the 
electron beam-plasma instability, a basic problem in kinetic plasma turbulence [3]. 
When a weak electron beam is injected into a thermal plasma, electrostatic modes at 
the plasma frequency (Langmuir modes) are destabilized and, subsequently, amplified 

B 

In both cases, the field intensity experiences a sudden growth, which is followed by a 
non linear saturation towards a non-equilibrium plateau. This initial violent relaxation 
is governed by the the Vlasov equation, a partial differential equation which represents 
the continuous counterpart of the discrete dynamics. Based on the Vlasov description, 
one can quantitatively predict the average behaviour of the system at saturation [6] and 
derive a reduced formulation to track the detailed time evolution of the main collective 
variables (TJ [H [9] . According to this latter scenario, a significant number of particles 
experience a bouncing motion in one of the (periodically repeated) potential wells, and 
form a clump that evolves as a single macro-particle localized in space. The remaining 
particles populate the surrounding halo, being almost uniformly distributed in phase 
space between two sharp boundaries. 

In real devices, however the finite extension of the interaction region, e.g. the size 
of the undulator, prevents the system to enter the deep saturated state and the initial 
sudden growth is the only regime that can be practically explored. It is therefore of 
general interest to mathematically address the study of the system dynamics for short 
times, aiming at providing closed analytical expressions that apply to a large class of 
initial condition. 

X Travelling Waves Tubes (TWTs) [3] have been extensively adopted to mimic beam-plasma 
interactions. The amplification process in a TWT can in fact be described in term of an analogous 
Hamiltonian setting. 
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When initializing the system around an equilibrium condition, one can 
straightforwardly linearise the equations of motion and obtain an estimate that is shown 
to accurately agree with the numerics. However, non-equilibrium initial conditions are 
also experimentally relevant [10] and result in a substantially different dynamics. In this 
paper, we shall focus on this latter case and derive perturbative solutions for the broad 
class of spatially non-homogeneous water bag initial profiles. 

The paper is organized as follows : In Section [2] we introduce the one-dimensional 
model of a FEL amplifier [5]. The continuum limit is also discussed and the Vlasov 
model presented. Section [3] is devoted to presenting key assumptions of the perturbative 
calculations. Closed expressions for the main macroscopic quantities are derived in 
Section [5] and compared with numerical simulations. Finally, in Section [6] we sum up 
and draw our conclusions. 



2. From the Hamiltonian model to the Vlasov equations: the case of the 
FEL 

In the following we shall focus on the case of a Free Electron Laser. However, as 
previously noticed, the model is more general, and can be regarded as a paradigmatic 
example of systems with wave-particles interactions. In this respect, the conclusions 
of our analysis will apply to other physical contexts where the same basic mechanism 
holds. 

Under the hypothesis of one-dimensional motion and monochromatic radiation, the 
steady state dynamics of a Single-Pass Free Electron Laser is described by the following 
set of equations: 



^ = V (1) 

dz Pj ' 1 j 

^1 = _ Ae Wi _ A * e -ie 3 ; (2 ) 

• < 3 > 



where A = A x + iA y represents the wave vector potential, 8j stands for the phase 
of the electron j with respect to the ponderomotive wave, while pj is its rescaled energy. 
All are adimensional quantities and the reader can refer to [5] for a detailed account on 
the derivation of the model and an exhaustive connection with the physical parameters 
of the machine. We shall here simply recall that z is the rescaled longitudinal coordinate, 
inside the undulator, which essentially plays the role of time. In the following, for the 
sake of simplicity, we shall replace it by t. 

The above system of equations (N being the number of electrons) can be derived 
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from the Hamiltonian [jj: 

n 2 ry n 

^ = E| + 2 V^E sin ^-^' ( 4 ) 

3=1 3=1 

where the intensity I and the phase of the wave are given by A = yTjN exp(— i<p). 
Here the canonically conjugated variables are (pj, 9j) for 1 < j < N and (/, 0). Besides 
the "energy" H, the total momentum P = . pj + / is also conserved. 

In the continuum limit Eqs.Q-Q are rigorously mapped into the following system 
of partial differential equations: 

— =p—- 2{A x cos6- A y sm0)—, (5) 



dA x 
dt 

dt 



d9 / dp f {6, p,t) cos 9, (6) 



oo 

oo 



d9 / dp f (6, p,t) sin6, (7) 



where f{0,p) represents the single particle distribution. The bunching coefficients 
bk are in turn defined as: 

del dpexp(-ik6)f(6,p,t) for k = 1,2,3.. (8) 

-7T J — OO 

and measure the degree of spatial packing of the particles. 

Assuming periodic boundary condition, the system ©-([7]) admits the following 
stationary solution: 

A x = A y = f = f (p) (9) 

Linearizing around equilibrium one can derive an approximate solution that holds 
for relatively short times. To this end, the following ansatz is put forward: 

f(9,p,t) = fo(p)+fx(6,p,t), A x (t)=X l (t) and A y (t) = Y l (t) .(10) 

Introducing in system ([5]) and using the notation i](p) = d p fo, we obtain at lowest 
order 



(dt+pde)h 



d6 I dp fi cos 9 — 

7T J — OO 

r p+oo 

d6 1 dp fi sin 9 + 



+ r?(Xi cos# - Vising = 
dX 1 



7T J — OC 



dt 
dY x 
~dT 








(11) 
(12) 

(13) 



§ In the following we set the detuning parameter to zero, thus assuming perfect resonance condition. 
The analysis can be however extended to the case where an energy mismatch has to be accounted for. 
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Such a linear system admits the following normal modes solution 



fi(9,p,t) = F x (p) S Q -^ + F*(p) e -^ 9 -^ 
X x (t) = X x e~ iuJt + X* e iuH 
Y x {t) = %Y X e~ iuit - iY* e iu>H . 



(14) 
(15) 
(16) 



where the symbols * stands for the complex conjugate and uo G C. By introducing 
solution (|T4 l) -f |T6l) into the linearized system one obtains, after some algebra, the 
following dispersion relation: 



which can be solved with respect to u, once the equilibrium initial condition fo(p) is 
specified. If a complex solution exists, the field grows exponentially, see Eqs. (IT4l-(IT6"l)). 
otherwise it oscillates indefinitely (see also [5]). 

It should be however stressed that the above treatment applies if the system is 
locally pertubed around the equilibrium initial condition (Q. For more general out-of- 
equilibrium settings, the linearization fails and other strategies need to be developed. 
In particular, we shall here discuss a perturbative approach aimed at characterizing the 
evolution of the system initialized in the so-called water-bag state, thus allowing for a 
spatial bunching (6^(0) 7^ 0) of the beam. This technique enables us to derive closed 
analytical expressions for the time evolution of the fundamental macroscopic observables 
that characterize the system dynamics.. 

3. Simplified water-bag approximations 

In the following we shall consider an initial water-bag profile: particles are confined in a 
finite portion of phase-space and there display a uniform distribution. This is a rather 
common choice already invoked in several studies [6] and often assumed to provide a 
simplified description of more realistic initial conditions. 

Our perturbative analysis is based on a simple assumption: we hypothesize that 
the initial particles' evolution changes the shape of the water-bag, while preserving the 
uniform distribution (homogeneous density) inside the stretched domain. This working 
ansatz is corroborated by direct numerical inspection and shown to hold approximately 
during the initial violent relaxation, until the saturated regime is eventualy attained. 
From there on, a dense cluster starts to develop and one has to resort to the so-called 
macro-particle scheme (see [8j [9]) to derive a reduced theoretical framework. 

According to the proposed formulation, and recalling that Liouville theorem holds, 
one can formally trace the water-bag evolution in term of its outer boundaries, once the 
initial density has been assigned. Figure [1] illustrates this concept pictorially, for the 
case of a rectangular water-bag domain which will be assumed in this study. Hereafter, 
the upper and lower boundaries are parametrized as P + (9,t) and P_(6*,t), a functional 
dependence which can be correctly invoked as long as the evolution stays single-stream, 




(17) 
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i.e. before a lateral flip occurs. As for the left (resp. right) lateral edge, we assume 
9 = Q-{t) (resp. 9 = + (t)), which in turn amounts to model it as a vertical, though 
dynamic, barrier, hence neglecting its inclination. These prescriptions translate in the 
following mathematical expression for the initial single-particle distribution function 



f(6, P , o) = f [Q(e - e_(t)) - n{e - e+(*))][ft(p - p_(0, t)) - q( p - p + (e, t)\, (is) 



where / labels the water-bag's density and Q represents the Heavyside 
function. Moreover, we shall limit the discussion to initially symmetric profiles which 
mathematically yields to: 



Notice that for + = — ©_ = n, one formally recovers the setting (jUJ) and an 
exponential growth is thus expected, provided A x = A y = 0. In the following we will 
concentrate on the case a < it. The normalization condition results in fo = l/2aAp. 

Particles and wave are in phase at t = 0, i.e. 0(0) = @+(0) — B_(0), a condition 
that results in the optimal growth of the field intensity initially set to Iq (either finite or 
zero). Finally, we shall assume that the P + and P profiles are accurately interpolated 
by two parabolas centered in 9 = 0, and consequently parametrized as: 



numerically via direct fit and further confirmed by a posteriori testing the predictive 
ability of the self-consistent theory here developed. The linear term in 9 does not appear 
in Eq. fTSO]) as the associated coefficient scales as t 3 , thus falling beyond the accuracy of 
our perturbative scheme. 

4. Solving the short-time system dynamics: A perturbative derivation 

Starting from this initial setting and recalling Eq.fjHJ), one straightforwardly obtains: 



f(o, P ,oy. 




(19) 




A x = fo(v+ - w_)(sin0 + - sin©_) = s a + 0(t) 



(21) 



where 



s a = sin(a)/a. Integrating yields : 
A x (t)=A° x + s a t + 0(t 2 ), 



(22) 



with A° x 




Iq. As for A y , combining Eqs. ([7j) and (1201) results in : 



A y = fo(v + — t>_)(cos 0+ — cos 0_) = 0(t), 



(23) 



since 0+(O) 



0_(O). This in turn implies 



A y {t) = 0(t 2 ). 



(24) 
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We now look at the evolution of a reference particle (here baptised test-particle) 
located at the boundaries of the waterbag. The motion of such a particle is governed 
by the following Hamilton equations : 

f 6 = p 

I p = -2(A,. cos - A y sm6) 1 ' 

Consider in particular a (virtual) particle of initial coordinates (9(0) = ±a,p(0) = 
0). During the time interval covered by our investigations, namely before the outer 
contour enters a multi- stream regime, it can be reasonably assumed that the particle 
evolves coherently with the associated boundary (see FigJT]), an observation which 
suggests identifying 6(t) = Q±(t). According to Eqs. (l22l - (l25l) . the particle's position 
obeys to : 



± = - 2(A X cos 6 ± - A y sin Q ± ) (26) 
~ - 2 A° x cos a + 0(t). 

The preceding equation can be integrated and results into 0±(t) = ±a — 
A° x cos (a)t 2 + 0(t 3 ). This latter expression is then re-inserted into the first of ( 1261) : using 
again relation (T22]) . after integration, one ends up with the more accurate expression : 

e±(t) = ±a- A° x cos (a)t 2 - \s a cos (a)t 3 + 0(t 4 ). (27) 

Following the same reasoning, we shall now consider particles initially positioned 
in correspondence of the upper, alternatively lower, boundary of the rectangular water- 
bag profile, (0(0) = 0,p(0) = ±Ap/2). These particles are also virtually linked to the 
boundary that they contribute to create and thus p(t) = P±(9(t),t). Recalling that the 
phases evolve as 9(t) = ±tAp/2 + 0(t 2 ), one can solve Eq.( l25l) to obtain the following 
expression for the conjugate momenta : 

P ± (9, t) = ±^ - 2A° x t - s a t 2 + 0(t 3 ) (28) 

Then, as that P± = u9 2 + v±, and since u(t) goes as 0(t), one can conclude that : 

v± = ±^ -2A° x t- s a t 2 + 0(t 3 ). (29) 

The conservation of the total momentum which, in the Vlasov picture, reads 
P = A 2 + Ay + jj f(9,p, t)p d9 dp, now takes the form : 

Jo = A\ + A\ + / (^(e 3 + - & 3 _)uAv + i(0 + - e_)tJA«), (30) 

where v = v + + w_ and Av = v + — w_ have been introduced. Moreover from 
Eqs. (J22J and ([24]) that : 

A 2 X + A 2 = 1° + 2A° x s a t + s 2 J 2 + 0(t 3 ). (31) 
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Expression (1317]) can be hence manipulated by making also use of expansions (1271) 
and fl29|) . Solving for u{t) results in : 

2a z a z a z 

Furthermore, the conservation of the energy reads: 

,2 



u{t) = —(s a -l)v + 0(t 3 ) = (l- Sa )A x t + —s a (l-s a )t 2 + O(t 3 ).(32) 



H = J J f(6,p,t)?j d9 dp+2 J j f(9,p,t)(A x sm6+A y cos 6)d6 dp, (33) 



which leads to : 

^ = _ e *_y Av + (9 3 _ e i )u7jAv (34) 

+ (6+ - Q^)^{Av 2 + 3v 2 )} + 2{A y A x - A x A y ). 

The l.h.t. is the initial energy of the system: since the initial water-bag is 
centered around zero, the field does not contribute to the energy at t = and the 
only residual component comes from particles' kinetic energy. Finally, replacing in the 
above expression each term - apart from A y - by its expansion (see equations (|27j) . fl29l) 
and (132]) ). we get for A y a first-order differential equation. Then, assuming the general 
form A y {t) = rjt n + ut n+l + 0(t n+2 ), one immediately realizes that n = 3, and, more 
precisely : 

A y (t) = ^(4 - 8s a + 9s 2 a )t 3 + ^(4 - 8s a + 9s 2 a )t 4 + 0(t 5 ). (35) 

Note that the expansions here derived can be used as a starting point to calculate 
higher order corrections, following a typical strategy often employed in perturbative 
analysis. In particular, from Eqs.Q, (127]) and (129]) . we get : 

A x {t) = A° x + s a t + 0{t A ). (36) 

It should be however stressed that this procedure cannot converge indefinitely, since 
the assumptions built into the model will eventually prove inaccurate and further effects 
will need to be properly incorporated (e.g. the leaning of the lateral boundaries). 



5. Predicting the macroscopic observables 

In the above paragraph we have developed a perturbative approach that ultimately 
enabled us to provide closed analytic expressions for the time evolution of the complex 
field A and the single-particle distribution function f(8,p,t). This novel insight 
allows us to condensate in compact formulae the time dependence of all fundamental 
macroscopic quantities, modified during the self-consistent amplification process. Few 
examples are discussed in the remaing part of this Section. 
The wave intensity I(t) follows trivially as: 

I(t) = \A\ 2 = A x (t) 2 + A y (t) 2 = I + 2y^T s a t + s 2 J 2 + 0(t A ). (37) 

For I = 0, the laser intensity scales quadratically with time, a result previosuly 
reported in [TT]. This finding agrees with direct numerical simulations based on 
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the N-body model (j3J and reported in figure [2J Clearly, our solution is limited to 
non- homogeneous (bunched) initial beam: If the phases of the particles are initially 
occupying the whole interval [— 7r,7r], an exponential instability develops, as predicted 
by the linear analysis of section [2j 

Consider now the case Iq ^ and define the gain G(t) = I(t)/I . One can hence 
recast equation fl37j) in the form : 




G(t) = l + 2-+ - +0 - . (38) 



where we have introduced the characteristic time T c = y/T^/s a . In principle, 
increasing T c amounts to slow down the growth: longer times (undulators) are thus 
required to attain a fixed gain level. Interestingly, a faster evolution is produced when 
increasing the initial particle bunching. A similar effect is obtained by reducing the 
intensity Io of the injected seed. In figure El we report the gain G as funtion of the 
rescaled time t/T c : symbols refer to numerical simulations based on Hamiltonian (j4j). 
The data are nicely interpolated by the universal profile fl38|) . By inverting equation 
( l38l) . one can estimate for time t* needed to the system to reach a fixed gain amount 
G*. A straightforward calculation leads to the following compact relation 



t* = y-^(VG*-l), (39) 
s a \ J 

which can be used as a first rough guideline for optimization and design purposes. 
As a second example, we consider the particles' energy dispersion here defined as : 



D = J J old dp f p 1 - (J J d6 dp f p)\ (40) 
From the above, after some algebra, it follows : 

D{t) = ^- + ^ll^ 



Ap 2 16 r2 ( Sa -l) 2 
12 + 5 si 




(41) 



For small values of Iq, the wave intensity increases as s 2 a t 2 and hence the particles' 
energy scatters as t 4 : 

An 2 1 

D{t) = ^ + -(44 - 8sl + As 2 a )t 4 + 0(t 5 ). (42) 

Again, the theory agrees well with direct simulations as reported in Figj4] (b). 
Finally, the bunching coefficients ([HD read : 

b k (t) = S ^^- + 0(t% for k = l,2.. (43) 
ka 

This in turn implies that there is no significant change in the bunching of the 
particles during the early stage of the interaction. 

As a side remark we shall notice that the higher order corrections for A y as predicted 
by equation (I3^|) correlate very well with the numerics, see Fig. 
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Figure 1. Left : Evolution of the bunch of particles in phase-space, at t = 
0, 0.25, 0.5, 0.75, 1, 1.25, 1.5 (from top right to bottom left). The boundaries 
P± are plotted in dark (horizontal boundaries), while the inner part of the waterbag 
stands in grey. At t ss 1.4, the bunch flips : this is clearly a limitation (in time) of 
our modelization of the system. Right : waterbag at t = 0.5; the dash-dotted lines 
correspond to the second degree polynomial fit of Eq. ([20]) . 




Figure 2. Time evolution of the laser intensity. Here Iq ~ 0. Particles are initially 
distributed in the interval [— f ; § ]• The solid line refers to the quadratic law predicted 
by Eq. (|37[) . The circles represents the numerical simulations based on Hamiltonian (TjQ) . 

6. Conclusion 

In this paper we developed a perturbative approach to characterize the short-time 
evolution of a Single-Pass Free Electron Laser. In particular we provide closed analytical 
expressions that elucidate the time dependence of the main macroscopic quantities, e.g. 
laser intensity, degree of bunching and energy dispersion. More specifically, particles 
are initially randomly distributed inside in- homogeneous (spatially bunched) water- 
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c 



Figure 3. The gain G is plotted as function of the rescaled time t/T c . The plain line 
refers to the theoretical prediction (|38p . the symbols to numerical simulations : the 
crosses correspond to Iq = 0.8N and a = tt/2, the circles to Iq = 0.8iVand a — 7r/4, 
the triangles for I = OAN and a — ir/2, while the stars stand for I = N and a = ir/2. 
We shall here notice that equation (f38|) is found to be accurate for initial bunching and 
intensity resp. smaller than tt/2 and 0.8N. For larger values, higher order corrections 
need to be incorporated into the model. 

0.4, , . . , 1 



0.3 
Q 0.2 




t/T 

c 

Figure 4. Diffusion D(t) vs time t. Symbols refers to the N-body simulation, while 
the dashed and solid lines represents the theoretical prediction (|4"Tj) : the former stands 
for the second order in t prediction, the latter for the third order. 

bag envelopes. The underlying Vlasov dynamics results in a progressive distorsion of 
the water-bag profile, which we here monitor by tracking the evolution of the lateral 
boundaries. The prediction of the theory are shown to agree with direct numerical 
simulations. Interestingly, we also derive a universal relation which allows one to 
calculate the laser intensity gain at a given undulator length. 

In conclusion, it is worth emphasising that we here address the study of a generic 
wave-particles interaction process. It can be thereofore expected that our conclusions 
will prove useful beyond the realm of FEL applications, and possibly translate to other 
fields where the complex interplay between particles and waves is known to be crucial. 
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0.8 
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Figure 5. A y (t) vs time t. Symbols refers to the N-body simulation, while the solid 
line represents the theoretical prediction ([55)1 . 
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